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r^ I momentum. We then turn to the magnetorotational instability, which thus far is the only 

mechanism that has been shown to initiate and sustain turbulence in disks. Finally, we 



2 : 1 Introduction 

It is an observational fact that at least part of the mass in protostellar disks is accreted onto 

.!^ ] the central objet. In a Keplerian accretion disk, the specific angular momentum increases with 

rS ' radius. Therefore, a particle can be accreted only if its angular momentum is removed or 

d • transferred to particles located at lager radii. Whether angular momentum is removed from 

or redistributed within the disk depends on whether the disk is subject to external or internal 

torques. Possible external torques can either be magnetic (when an outflow is present) or tidal 

(in binary systems), whereas possible internal torques can either be gravitational (massive 

disks) or turbulent. These mechanisms have been discussed by Papaloizou & Lin (1995). 

During the early stages of disk evolution, when the disk is still embedded (class 0/1 object) 
and has a significant mass compared to the central star, there may exist strong disk winds and 
bipolar outflows (e.g. Reipurth et al. 1997) with associated magnetic fields. During this stage, 
a hydromagnetic disk wind may be an important means of angular momentum removal for the 
system (Blandford & Payne 1982). Because of the action of magnetic torques, material ejected 
from the disk is able to carry away significantly more angular momentum than it originated 
with in the disk. Therefore, even a modest ejection rate can lead to a significant accretion 
rate through the disk. However, observations indicate that outflows may exist only in the 
early stages of disk evolution, so that this mechanism cannot account for angular momentum 
transport during the whole life of accretion disks. In addition, it may affect only the inner 
parts. 



When the mass of the disk is significant compared to that of the star, gravitational instabihties 
may develop, leading to outward angular momentum transport (Papaloizou & Savonije 1991; 
Heemskerk et al. 1992; Laughlin & Bodenheimer 1994; Pickett et al. 1998) that results in 
additional mass growth of the central star. This redistribution of mass may occur on the 
dynamical timescale (a few orbits) of the outer part of the disk and thus may be quite rapid: 
on the order of 10^ yr for a disk radius of 500 AU. The parameter governing the importance 
of disk self-gravity is the Toomre parameter, Q ~ M^H/{M(ir)^ with M* being the central 
mass, M(i being the disk mass contained within radius r and H being the disk semi-thickness. 
Typically H/r ~ 0.1 (Stapelfeldt et al. 1998) such that the condition for the importance of 
self-gravity, Q ~ 1, gives M^ ~ O.IM,.. During the period of global gravitational instability, it 
is reasonable to suppose that the disk mass is quickly redistributed and reduced by accretion 
onto the central object such that the effects of self-gravity become negligible. 

If the disk surrounds a star which is in a pre-main sequence binary system, tidal torques 
transport angular momentum outward, from the disk rotation to the orbital motion of the 
binary. However, although tidal effects are important for truncating protostellar disks and for 
determining their size, it is unlikely that tidally-induced angular momentum transport plays a 
dominant role in the evolution of protostellar disks (see Terquem 2001 and references therein). 
In a non self-gravitating disk, the amount of transport provided by tidal waves is probably 
too small to account for the lifetime of protostellar disks. In addition, tidal effects tend to be 
localized in the disk outer regions. 

When the disk mass is such that self-gravity can be ignored and the jet activity has significantly 
decreased, turbulent torques may become the most important way of redistributing angular 
momentum in the disk. Historically, the first angular momentum transport mechanism to be 
considered was through the action of viscosity (von Weizsacker 1948). However, in order to 
result in evolution on astronomically interesting timescales, it is necessary to suppose that an 
anomalously large viscosity is produced through the action of some sort of turbulence. 

In this paper we present the theory of turbulent accretion disks. We first calculate the molecular 
viscosity in a disk, and show this is far too small to account for the evolutionary timescale 
of disks. We then describe how turbulence may result in enhanced transport of (angular) 
momentum. About ten years ago, a linear magnetic instability was identified (Balbus & Hawley 
1991) that results in MHD turbulence. We present this instability. We then describe both the 
basis and the structure of a disk models, which were introduced by Shakura & Sunyaev (1973). 



2 Fundamental of hydrodynamics 

2.1 Fluid mechanics vs. kinetic theory 

The fundamental equations of fluids can be derived by considering them as either a collection 
of particles (kinetic theory) or as a smooth continuum. This latter approach is justified when 
the mean free path A of the particles is very small compared to the macroscopic length scale L 
of interest in the fluid. This condition is not met in stellar systems (galaxies or star clusters), 
where A 3> L, nor in planetary rings, where A ~ L, but it is fulfilled in disks. The conservation 
equations are obtained more straightforwardly when considering the fluid as a continuum, which 
is the approach we will use here. 

This does not lead to a precise expression for the transport coefficients, however, in contrast 
to kinetic theory. Transport of energy, mass and momentum occurs in a gas which is out of 



equilibrium (i.e. in which the distribution function is not a Maxwell-Boltzmann distribution) 
through molecular collisions. Most of the time, the departure from equilibrium is tiny, so that 
the distribution function is nearly maxwellian. Within the context of the kinetic theory, in 
which molecular collisions are explicitly calculated, the Chapman-Enskog procedure gives the 
transport coefficients by considering small variations of the distribution function around the 
Maxwell-Boltzmann distribution. Such a calculation is not possible when fluids are viewed as 
continua, as in this case molecular collisions are not explicitly calculated. It is possible however 
to get a phenomenological expression for the transport coefficients in this context, as we shall 
see below. 

For more details on the kinetic theory of gases, the reader is referred to, e.g., Huang (1987) and 
Shu (1992). For the description of fluids as continua, which is developed below, a more complete 
presentation can be found, e.g., in Acheson (1990), Batchelor (1967), Landau & Lifchitz (1987) 
and Shu (1992). 



2.2 Conservation laws 

2.2.1 Mass conservation 

If there is no sink or source of matter, the time rate of mass variation within a flxed volume 
V is equal to minus the flux of mass (mass advected by the flow) through the surface S which 
encloses the volume: 

where p{r,t) and v(r,t) denote the density and velocity vector of the fluid at location r and 
time t, and n is the unit normal to the surface oriented outward. 

Since this has to be satisfled for any volume V, and using the divergence theorem, we get: 

|^ + V-(pv)=0. (2) 

The mass conservation equation derived above is also know as the continuity equation. It states 
that mass variation in a flxed volume is due to a continuous flow of mass through the surfaces 
which enclose the volume, and not to 'jumps' from one location to another remote one. 

Equation (|^) has been derived by considering the variation of mass of a volume flxed in space. 
There are two contributions to the time rate of mass change in this volume: pV • v, which is 
due to the change of the volume itself, and v ■ Vp, which is due to the mass flowing through 
the volume. 

We can also consider the time rate of mass change of a volume of fluid that we follow along its 
trajectory. For that we deflne the Lagrangian derivative d/dt, which is the time rate of change 
in a frame moving with a particular fluid element along its trajectory. The density p depends 
both on r and t and, if r is the location of the moving fluid element, it depends also on t. Then 

dp ^ dp df^dx^ ^ df^ 

dt dt^dx,dt at ^ ^' ^^ 



were the Xi {i = 1,2, 3) denote the coordinates and we adopt the standard Einstein notation of 
summation over repeated indices. 

Equation (|^) can therefore be rewritten as: 

|+pV.v = 0. (4) 

This expresses the fact that the mass change in a volume moving with the fluid is due only to 
the change in the volume. 

2.2.2 Incompressibility 

A fluid is said to be incompressible if a volume element is neither compressed nor dilated when 
moving with the flow. Let us consider a fixed volume V of the fluid bounded by the surface S. 
The volume of fluid which leaves V through an infinitesimal surface dS per unit time is v ■ ndS, 
where n is the unit normal to the surface oriented outward. Therefore 

dV r r 

^ = VndS= / (V • v) rf\/, (5) 

dt Js Jv 

where we have used the divergence theorem to get the final expression. The condition for the 
fluid to be incompressible is therefore: 

V ■ V = 0, (6) 

or, equivalently, using the mass conservation equation (^): 

This expresses the fact that, in an incompressible fluid, the density of a fluid element is constant 
along its trajectory. 

2.2.3 Equation of motion and momentum conservation 

The product of the acceleration of a fluid element moving with the flow, dv /dt, by the mass 
density p is equal to the net force per unit volume exerted onto this element. This has contri- 
bution from external forces (e.g., gravitational, magnetic.) and/or from momentum transport 
within the flow due to molecular collisions (pressure and viscous forces). The local form of the 
equation of motion is then: 

where P is the pressure, 'Fvisc is the viscous force per unit volume and F is any other force per 
unit volume which does not arise from pressure and viscosity. When "Fvisc = 0, this equation is 
known as Euler equation. 

Using equation (EI), we can rewrite the equation of motion under the form: 



- (pv) = -V ■ (/)vv) - VP + F,,,, + F, (9) 

where, according to our notation, the i component of the vector V ■ (pvv) is V ■ {pviv). 

This is the momentum conservation equation, which states that the time rate of momentum 
variation within a fixed volume is equal to minus the momentum which is advected by the flow 
through the surface which encloses the volume plus the net force exerted onto the volume. 



2.3 Viscous forces: Molecular transport of momentum 

The viscous force in a fluid is due to an irreversible transport of momentum from regions where 
the velocity is higher to regions where it is lower. A variational calculation shows that uniform 
rotation is an extremum for the energy. Sometimes it is a minimum, in which case dissipation 
of energy at flxed (angular) momentum results in uniform rotation (e.g. synchronous binaries). 
This tends to be the case when pressure is dominant against gravity. Sometimes the extremum 
is a maximum, in which case all the matter tends to accumulate at the center, with all the 
angular momentum at inflnity. This tends to be the case when rotation is dominant against 
gravity, and this is what happens in disks (for further details, see Lynden-Bell & Kalnajs 1972; 
Lynden-Bell & Pringle 1974). The dissipation of energy results from molecular collisions, or 
friction. 

We denote by —Uij the flux of the i component of the momentum in the j direction. The tensor 
[a] is called the viscous stress tensor. The i component of the viscous force per unit volume is 
then: 

F =^ nm 

J^ vise A r% • V-*-^/ 

dx, 

Friction occurs only when different parts of the fluid have different velocities. Therefore aij 
should depend on the velocity gradients. If the velocity varies on a scale large compared to 
the mean free path, i.e. to the scale over which molecular transport arises, one can suppose 
that aij depends only on the flrst derivatives of the velocity with respect to the coordinates. 
Furthermore, we suppose that the dependence is linear, i.e. we limit ourselves to Newtonian 
fluids. The most general form of aij is then: 

,,«|li + ^|i + B*,f^, (11) 

OXj OXi OXk 

where A and B are constants to be determined, and 6ij is the Kronecker symbol. If the flow 
is uniformly rotating in the {xy)-plane for instance, we must have a^y = 0. Since v^ cc y and 
Vy oc —X in that case, that implies A = 1. Therefore the tensor [a] is symmetrical. We note 
that: 

Tr [a] = (2 + 35) V ■ v, 

which shows that Tr [a] is a measure of the volume change of a fluid element. It is an experi- 
mental fact that the stresses which change the volume of a fluid element give different viscous 
forces than the stresses that preserve the volume. Therefore we rewrite aij under the form: 



where the term in bracket on the right hand side is trace free, i.e. does not modify the volume 
of a fluid element. 

The shear and bulk viscosities, that we denote r] and ( respectively, are then experimentally 
defined as: 



a., = r^ -^ + ^ - -<5., V ■ V + (S.^V ■ v. (12) 




The bulk viscosity is associated with internal degrees of freedom of the molecules in the fluid. 
It becomes negligible if the equipartition between these different degrees of freedom is reached 
over a timescale shorter than the timescale between two collisions. Furthermore, for a perfect 
monoatomic gas it can be shown that C = (Huang 1987). Therefore, from now on we shall 
ignore it. 

With this expression for aij, we can rewrite the i component of the momentum equation (^ 
under the form: 

{pViVj + Tij) = Fi, (13) 



dt dxj 



where 



T,, = PS,, - a.,. (14) 

Equation (p!3D makes is clear that the flux of the i component of the momentum in the j 
direction, fnJiVj + Tij, has contributions from both advection and molecular transport (pressure 
and viscous forces). 

For an incompressible fluid in which rj is constant, the corresponding equation of motion, derived 
from equation (^ with the above expression of the stress tensor, is called the Navier-Stokes 
equation. 

We note that since friction converts mechanical energy into heat, the rate of change of the 
kinetic energy, dE^/dt, must be negative if only viscous forces are acting. That can be shown 
from equation (|T3|) to be equivalent to the condition r] > (e.g.. Landau & Lifchitz 1987, 
§ 11.16). 

2.4 Expression of the shear viscosity 



The expression ([l^ for the stress tensor need not be exact. It has been derived phenomeno- 
logically assuming that a^j depends only on a linear combination of the first derivatives of the 
velocity with respect to the coordinates. However, the kinetic theory applied to dilute gases 
leads to the same expression for a^j. Here we are going to derive an approximate expression for 
the shear viscosity t] by using the kinetic theory of gases, and we will point out that it arises 
from correlations between the particle velocities. 
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2.4.1 Kinetic theory of gases 

For simplicity, let us assume that the velocity in cartesian coordinates is v = Vx{y)e.xi where 
Bx is a unit vector along the x-axis. The expression ([T2|) for the stress tensor then gives: 

dvx , . 

(Jxy = CTyx = V-^-- (15) 

dy 

This relation had already been proposed by Newton and Hooke, but it is Maxwell who gave 
the derivation of the viscosity r] that we develop now. 

In average, a molecule has a collision with another molecule after it travels through a distance A, 
which is the mean free path of the particles. We suppose that after the collision, the momentum 
of the molecule is the same as that of its new environment. 

Let us consider the momentum which is transported during the time St across a surface element 
6S perpendicular to the y-axis and with ordinate y. There are nuStSS/6 particles crossing 6S 
from above during 6t, where n is the number density of particles, u is their random (thermal) 
velocity relative to the mean flow and the factor 6 comes about because there are three possible 
directions for the particles, each with two orientations. Each of these particles travel through 
A before it suffers a collision below 6S, which results in its momentum varying by: 

m [vx{y) - Vx{y + A)] ~ -mX—^, 

to first order in X/L, where L is the scale of variation of the velocity. In other words, each 
particle carries below SS the excess of momentum mXdvx/dy. Here m is the mass of a particle. 
On the other hand, each particle traveling upward carries above SS the deficit of momentum 
—mXdvx/dy. Therefore, the net x-component of the momentum which is carried downward 
during 6t by the particles crossing dS is: 

6'^Px = 2 I -nu6t5S ] mX—^ = -nmuX—r^6S6t. (16) 

^ V6 J \ dt J 3 dy ^ ^ 

Let us now consider a box with horizontal faces at y and y + Sy and surface area 6S. The 
exchange of particles across the upper face during the time 6t results in the momentum S'^px{y+ 
6y) being added to the volume, whereas the exchange across the lower surface results in the 
momentum 6'^px{y) being removed from the volume. Therefore, the time rate of change of the 
momentum content of the box is: 



^S'^Pxiy + Sy) - S'^pxiy)] ^ j- i -nmuX—^ J 5S5y (17) 



5t 



to first order in X/L. This is also Fmsc,x SS6y, where Fmsc,x is the viscous force per unit volume. 
Therefore: 

Fyiscx — -r -nmuX—r^ . (18) 

' dy\3dyj ^ ' 

Since F^sc^x = daxy/dy and using equation (|T3p, this leads to: 
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V - -^PuX, (19) 

where p = nm is the mass density of the particles. We also define the kinematic viscosity 
V = T]/p, i.e. 

u ~ -uX. (20) 

The thermal velocity u is on the order of the sound speed. 

We denote by TZ the Reynolds number, which measures the relative strength of the inertial and 
viscous forces. The inertial force is ~ pv'^/L, where v is the flow velocity, and the viscous force 
is ~ puv/L"^ (from eq. [|T0| and [Q). Therefore, 



7^ = Lv/u. (21) 



2.4.2 Velocity correlation 



Here again, for simplicity, we consider the case where the flow velocity is along the x-axis and 
the motion is in the (x|/)-plane. The components of the instantaneous velocity of a particle 
in the fluid are {vx + Ux,Uy), where Ux and Uy are the components of the random (thermal) 
velocity relative to the mean flow. We have < Ux >=< Uy >= 0, where the brackets denote a 
time average. The flux of the x-component of the momentum along the y-direction is: 

nm {vx + Ux) Uy. 
Averaged over a large number of particles, or, equivalently, over time, this gives: 

P {UxUy) , 

since < VxUy >= Vx < Uy >= 0. By construction, this quantity is a component of the stress 
tensor. Therefore, 



a 



xy 



-P {UxUy) , (22) 



where the minus sign comes about because of the deflnition of the viscous stress tensor (see 
§ p.3| ). For a Maxwell-Boltzmann distribution function (or any xy symmetric function), < 
UxUy >= and there is no transport. Equations (p^D and (|^) lead to: 

{UxUy} = -^-j-- (23) 

Note that in general the time average of the fluctuation velocity may not be zero. In an accretion 
disk for instance, there is a net radial drift of mass, i.e. the mean value of the radial component 
of the fluctuation velocity is flnite. Here, ii < Uy > were non zero for instance, the flux of 
the x-component of the momentum along the y-direction would have the extra contribution 
pVx < Uy >. This represents the transport of mean momentum by the fluctuations. Formally 
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however, what we call the stress tensor would still be given by equation (p2D, as this is the only 
contribution from friction between particles. 

We denote by C the correlation coefficient between the velocities Ux and Uy-. 

{UxUy) 



c = 



V? 



where u is the characteristic velocity of the particles in the fluid relative to the mean flow. We 
see from ( ^2]) that C gives a measure of the stress tensor. If Cg is the sound speed, then m ~ c^. 
Therefore 

v dvx Xvx A 



cl dy L Cs L 



where we have used equation (|20|) and dv^/dy ~ v^/L. Here A4 = v^/cs is the Mach number. 
Using equations ( PD| ) and (0), we can write TZ = M.L/X, and therefore C = M'^/TZ. 

In an accretion disk, A4 ~ 10-20 and the Reynolds number corresponding to the molecular 
viscosity is larger than 10^^, so that C ~ 10^^^ ^ 1!! 

This is a consequence of the small value of the ratio of the mean free path to the scale of the 
mean flow (~ 10~^^) and it means the state of the gas is not affected at all by the molecular 
transport of angular momentum. In other words, the flow and the molecular transport are 
completely decoupled. 



3 Angular momentum transport by turbulence 

The realization that molecular transport of angular momentum is so inefficient led the theorists 
to look for another mechanism of transport in accretion disks. Because Reynolds numbers are so 
high, it was thought that probably accretion disks would be subject to the same hydrodynamical 
nonlinear instabilities that shear flows experience in laboratory. The resulting turbulence would 
then transport angular momentum efficiently. Although today much doubt has been cast on 
hydrodynamical instabilities in disks, turbulence is still a strong candidate for transport since 
it has been shown relatively recently that a linear magnetohydrodynamical instability can 
develop in disks (see below). Therefore, we turn now to turbulent transport, and contrast it 
with molecular transport. Much of this section is based on Tennekes & Lumley (1972), which 
the reader is referred to for more details (see also Balbus & Hawley 1998). 

We shall restrict ourselves here to the study of incompressible flows, as this simplifles the 
discussion. For our argument, it is sufficient to take into account only pressure and viscous 
forces, but in principle any other (external or inertial) force could be added. The equations 
describing the fluid are: 

^+^.^ = 1^^.. (24) 

dt ^ dxj pdxj *'^' 

and 



where the tilde symbol above a variable means we consider the instantaneous value of the 
variable at the location Xi and time t. Here 

dij = -p6ij + r]Sij, 

with 

dvi dvj 

^^ dxj dxi 

and p is the pressure. We suppose 77 is a constant. 

We now use the so-called Reynolds decomposition, in which an instantaneous value is written 
as the sum of a mean value (denoted by a capital letter) plus a fluctuation (denoted by a small 
letter): 

Vi = Vi + Vi. 

This fluctuation is characteristic of the turbulence. This decomposition is meaningful only if 
the timescale on which the fluctuations vary and the evolution timescale of the flow are well 
separated. The mean values are then taken over a timescale large compared to the turbulence 
timescale but short compared to that of the flow evolution. As here we are not interested in 
the long term evolution of the flow, we neglect the derivative with respect to time of the mean 
values (i.e. we consider a quasi-steady state). To simplify the discussion, we suppose that the 
average of Vi over time is zero: < Vi >= 0. 

Equation (|25| ) averaged over time then leads to dVi/dxi = 0, i.e. the mean flow is incompress- 
ible. Equation ( pSj) thus implies dvt/dxi = 0, i.e. the fluctuations are also incompressible. 

Using 

with (aij) = 0, the equation of motion averaged over time gives: 

^ -y^v^ + -ir i^j^^) = -4r^^v (26) 



dxj dxj pdxj 

where we have used dVi/dt = and the incompressibility of the mean flow and the fluctuations. 
The term < VjVi > represents the averaged transport of the fluctuations of the momentum by 
the fluctuations of the velocities. This is the turbulent transport. It is non zero only if the 
turbulent velocities in the different directions are correlated. It is an experimental fact that this 
is in general the case for shear turbulence (see below). Equation ( p6D shows that momentum is 
transferred between the fluctuations and the mean flow through the term d < VjVi > /dxj. 

We can rewrite equation (p6|) 

a 1 A 

(27) 

where Tjj = —p (vjVi) is called the Reynolds stress tensor, or turbulent stress. This is the 
contribution from the fluctuations to the averaged total stress tensor Sjj + Tij. Note that if 
we had not supposed < Vi >= 0, we would also have had terms like Vj < Vi >, representing 
transport (or advection) of mean momentum by the fluctuations. This is what Lynden-Bell & 
Kalnajs (1972) call lorry-transport, because it is a direct "shipment" by the equilibrium flow. 
Formally, this term is not part of what we call the turbulent stress however. 
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under the form: 










al/'^-- 


1 
P 


d 

dxj 


(S., 


+ 


r^J 



As we have no expression for the components of Tjj (six of which are independent), the problem 
has more unknowns than equations. This is the well-known closure problem of turbulence. 
Since Tjj appears in the same way as Hij in equation (pTj), it is tempting to express Tjj by 
analogy with Sjj, which depends on the molecular motion as we have seen in the previous 
section. This is the basis for the mixing length theory, in which Tij is written exactly like the 
tensor deriving from molecular motions using a so-called turbulent viscosity ut- By analogy 
with expression ( pOD for the molecular viscosity, it is supposed that z/^ ~ fyA, where vt is a 
characteristic velocity of the turbulent eddies and A is the so-called mixing length, which is 
the "mean free path" of the eddies, i.e. the distance they travel through before they mix with 
their environment. 

The same analogy is used in accretion disk theory through the a model that we shall describe 
in details below. 

It is important to note that the basis for this analogy is very weak. For a thorough discussion 
of the differences between molecular and turbulent transport, we refer the reader to Tennekes 
& Lumley (1972). Here we recall the main points. 

There is a fundamental difference between molecular and turbulent transport, in that turbulence 
is part of the flow whereas molecules are part of the fluid. We saw above that molecular 
transport is decoupled from the state of the flow. This is not true for turbulent transport, 
which depends completely on the flow for its existence. Turbulence can be sustained only if 
the eddies are able to tap the energy of the mean flow. Observations of laboratory shear flows 
for instance indicate that the eddies which are most efficient in tapping the energy of the mean 
flow are also most efficient in preserving a good correlation between the different components 
of the turbulent velocities, thus giving a large stress tensor. This is very different than for the 
molecules, whose velocities do not depend on non local properties of the mean flow. 

And as a matter of fact, transport of momentum is observed to be a characteristic of turbulent 
shear flows. In the air for instance, the correlation coefficient between the velocities is on the 
order of 10~^ for the molecules and 0.4 for turbulent motions. This is a characteristic of most 
shear flows. 

In addition, we saw when we derived the expression of the viscosity (|T9|) that this requires 
\ <^ L. For a turbulent flow, that would mean that the scale of the turbulent eddies is very 
small compared to the scale over which the characteristics of the flow, like the mean velocity, 
vary. This is not necessarily satisfied in an accretion disk, where both may be on the order of 
the disk semi-thickness. 

The conclusion is that shear turbulence is indeed a very efficient way of transporting momentum, 
but great care should be taken when expressing the resulting stress tensor by analogy with 
molecular transport. 

In particular, while representing gross transport by an effective viscosity can often be useful, 
doing a detailed stability analysis on a viscous fluid model for turbulent flow is generally not 
self-consistent, and can be very misleading (e.g. Hawley, Balbus, & Stone 2001). 

Note that the above discussion applies to shear turbulence only, i.e. flows where the fluctua- 
tions get their energy from the mean velocity gradients. Transport of momentum may not be 
nearly as efficient when the energy source for the fluctuations is a gradient of temperature or 
magnetic fleld for instance. As a matter of fact, there are strong indications that the transport 
of momentum associated with thermal convection is orders of magnitude weaker than that 
associated with shear turbulence (Balbus & Hawley 1998 and references therein). 
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Although turbulence has been considered as a way of transporting angular momentum in ac- 
cretion disks for more than fifty years, it is only relatively recently that an instability able 
to extract the energy of the shear and put it in the fluctuations has been elucidated. This 
instability, which requires a magnetic field, is described in the next section. 



4 Magnetohydrodynamic instabilities 

We consider now the stability of a disk containing a magnetic field. Balbus & Hawley (1991, 
1998 and references therein) have shown that in under very general conditions a magnetized 
disk is subject to a linear instability which nonlinear development was subsequently found to 
result in turbulence (Hawley, Gammie & Balbus 1995; Brandenburg et al. 1995). Following the 
presentation by Balbus & Hawley (1998), we first briefly describe the waves which propagate 
in a fluid containing a magnetic field (more details can be found in, e.g., the textbooks by 
Jackson 1975, Shu 1992 or Sturrock 1994), and we then go on to describe the linear instability. 



4.1 Waves propagating in a magnetized fluid 

Here we take into account only pressure and magnetic forces, as these are the only important 
forces for our argument. The governing equations are the induction equation: 

— = Vx(vxB), (28) 

which ensures that V ■ B = 0, the equation of continuity: 

f^+V- (pv) = 0, (29) 

and the equation of motion: 

dv 11 

— + (v-V)v = — VP + (VxB)xB, (30) 

dt p Hop 

where the last term on the right hand side is the Lorentz force per unit volume in SI units {po 
is the permeability of vacuum). Note that it can also be written: 

— (VxB)xB = -V|— 1 + I — • V)B (31) 

which shows that the magnetic force if the sum of a magnetic pressure and a magnetic tension 
(first term and second term on the right hand side, respectively). 

We first consider a nonrotating fluid and we look for small perturbations to an equilibrium state 
characterized by the quantities po, vq = 0, Pq and Bq. We denote the Eulerian perturbations 
by a prime, so that p = po + p' etc, and we suppose they vary on a scale much smaller than the 



scale of variation of the equilibrium quantities. The linearized equations (28), (29) and (30) 
are then: 
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-^ = Vx(v'xB„), 

dt po /ioPo 



;VxB') xBo. 



(32) 
(33) 
(34) 



Since the coefficients of the above equations are constant, we perform a Fourier transform for 
the variables r and t and look for solutions under the form exp[i {ut — k ■ r)] with u and |k| 
real constants. The above equations can then be rewritten by replacing d/dt by iu and V by 
-zk. 

To close the system of equations, we need an equation of state for the fluid. We consider 
adiabatic perturbations (i.e. we suppose there is no dissipation of heat), so that: 



P' p' 
P p 



(35) 



where 7 is the ratio of specific heats (7 = 5/3 for a perfect gas). 

The system of equations (p^)-(|35|) gives a dispersion relation of order six. It can be simplified 
by noting that the solutions associated with motions perpendicular to the (k, Bq) plane and 
those associated with motions in this plane can be decoupled. 

The ffist type of solutions is an Alfven wave, and its dispersion relation is: 



to = k Vj^ cos ■?/', 



(36) 



where ■?/' is the angle between k and Bq and 



va 



Br 



V^oPo 



is the Alfven speed. These are transverse waves (with no motion in the (k, Bq) plane) which 
propagate along the field lines. They do not involve any compression across the field lines and 
their restoring force is the magnetic tension. These waves are analogous to the waves which 
propagate along a stretched string. Within the factor cos-?/', we see that the Alfven speed is 
the square root of the magnetic tension divided by a mass density, just like the phase speed in 
the case of a string. 

The second type of solutions is associated with motions only in the (k, Bq) plane. Its dispersion 
relation is: 



u 




+ c ± 



{v\ + cl)' 



— 4:v\cl cos^ ip 



nl/2' 



(37) 



where c^ = •yPo/po is the adiabatic sound speed. The upper sign, which gives a larger phase 
speed uj/k, corresponds to the fast MHD wave, whereas the lower sign corresponds to the slow 
MHD wave. The fast mode is also referred to as magneto-acoustic wave. If either c^ <C va, 
va <^ Cs or cos-?/" <^ 1, then 

col = k' {v\ + c^) 
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and 



2 _ k^v\cj. cos^i) 

^- ^ 2~: 2 ' 

n + ct 



where the subscripts + and - denote the fast and slow waves, respectively. For the fast mode, 
the magnetic and thermal pressures act in phase. If Vj^ ^ Cg, the slow mode is an acoustic 
wave propagating along the field lines, whereas if va ^ Cg it degenerates into an Alfven mode 
in its dispersion properties (the eigenvector is distinct from that of the Alfven mode however, 
as the motions are not in the same plane). For the fast mode it is the opposite. 

In the absence of rotation, these modes are stable. However, Balbus & Hawley have shown that 
when rotation is introduced, the slow mode can become unstable, and this what we describe 
now. 

4.2 Linear magnet orot at ional instability 

We consider the simplest system in which the instability can develop. This is an axisymmetric 
disk with a vertical uniform magnetic field. For the original presentation, which includes the 
case of a radial field, see Balbus & Hawley (1991), and for the stability of a toroidal field see 
Balbus & Hawley (1992b), Fogglizo & Tagger (1995), Terquem & Papaloizou (1996) and Ogilvie 
& Pringle (1996). 

Since it is the slow mode which is destabilized, one can consider an incompressible fiuid. That 
amounts to taking the limit c^ — > cxd in the above equations, so that the fast mode disappears 
(its frequency becomes infinite). 

The system of equations describing the fiuid is then 

— = Vx(vxB), (38) 

V ■ V = 0, (39) 

(9v 11 

— + (v-V)v = — VP + (VxB)xB-V^, (40) 

dt p flop 

where \E' is the gravitational potential due to the central star. 

We use the cylindrical coordinates (r, 0, z) and we denote the unit vectors in the three directions 
by e,,, e,^ and e^. The components of the vector (v • V) v which appears in the equation of 
motion are: 



Or r d(f) dz r ' 

'"^^T + ~id + ^^ir^ + — ^' 

or r 0(p oz r 

dvz , ^^ , ^ 
dr r d(j) dz 

The equilibrium quantities, that we suppose uniform, are Pq, Pq, and Bq = -Bq^z- The velocity 
is Vq = rVL{r)e^. We now suppose that this equilibrium is slightly perturbed, so that p = Po + p' 
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etc, where the prime denotes the Eulerian perturbations, and we look for solutions proportional 
to exp[i{ujt — kzZ — krv)]. Here we consider axisymmetric perturbations, but more general 
solutions can be obtained (see the above mentioned papers). 

The linearized system of equations is then: 



iuBl = -ihBvl, (41) 

iuB'^ = -ik,Bv'^-ik,rQB'^+l-^-ikrrQ\B'^, (42) 

iuB' = v'^ + ikrBv'^, (43) 

r 

- - ikr] v' - ikv' = 0, (44) 

r / 

. ; ^^ , ikrP' ik.BBl iKBB'^ . ^. 

lujv'^ - 2VLv'^ = — ^ + — ^, (45) 

P Pop Pop 

..., + (^^^ + fij .. = --^^-^, (46) 

^u:v'^ = ^, (47) 

where, for simplicity, we have dropped the subscript '0' for the equilibrium quantities. We now 
consider large vertical wavenumbers, such that \kj^ 3> \kr\ and |fcz| ^ 1/r. Then, the above 
system of equations leads to the following dispersion relation: 

'^'-^(2*\U^ + 4tf)+.M(.M + ^)^0. (48) 

There is instability if uj has a negative imaginary part (the perturbation then grows exponen- 
tially with time). Equation (^) is a quartic for uo'^ which solutions are real. Therefore there is 
instability if uo'^ < 0, which requires: 

ev\ < -^. (49) 

This criterion has a very simple physical explanation. It states that there is an instability when 
the magnetic tension that acts on a segment of a field line is smaller than the net tidal force 
(i.e. centrifugal force minus gravitational force) acting on it. 

For a given equilibrium field B, and therefore a given Alfven speed va-, there will always be 
a wavenumber k such that this inequality is satisfied provided the right hand side is positive. 
Therefore, all the disks with 

are unstable, and this is the criterion for instability. 

The heart of the instability resides in the fact that a perturbed fluid element tends to conserve its 
angular velocity when a magnetic field is present. This is to be contrasted with a non magnetized 
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disk, in which a perturbed element tends to conserve its specific angular momentum. When 
displaced inward therefore, it has too much angular momentum for its new location (as the 
angular momentum increases outward in an accretion disk), and it moves back to its initial 
unperturbed position. When a magnetic field is present, the magnetic tension along the field 
line tends to enforce isorotation of the elements to which it is connected. A fiuid element 
displaced inward has therefore a lower angular velocity than the elements at its new location, 
and thus not enough angular momentum for its new position. As a result it sinks further in. On 
the opposite, a fiuid element connected to the same field line and displaced outward will tend 
to move further out. Angular momentum is transferred via magnetic torques from the inner 
fiuid element to the outer fiuid element. Note that the source of free energy for the instability is 
not in the magnetic field, but in the disk differential rotation. The magnetic field just provides 
a path to extract the energy. 

From equation (|48D, we can write the negative values of uJ^ as a function of k'^v\, and show 
that the maximum growth rate is: 



kc;n 



dVL 



d\nr 



(51) 



For a Keplerian disk. 



k^n 



?n 



(52) 



and is attained for 



kv. 



Q. 



(53) 



This holds even if the field has radial and azimuthal components and if the perturbed quantities 
are allowed to vary with r and (j) provided we then replace kvA by k ■ v^. Note that the non- 
axisymmetric case is more subtle though, as in that case plane waves cannot be sustained, 
being sheared out by the differential rotation. If we write k = kr^r + {'m/r)e 
being the azimuthal wavenumber, then kr is time dependent and 



+ kzGz, with m 



krii) = ko — mt 



dft 
dr 



(Goldreich & Lynden-Bell 1965), which means that a disturbance always becomes trailing in a 
disk where the angular velocity decreases outward. If kr is initially large and positive (leading 
disturbance), then the mode is stable as indicated by (|49|). But as time goes on, k^ decreases, 
so that k enters a region of instability. As kr becomes negative and keeps decreasing though, 
the mode becomes stable again. Formally, the instability is therefore not purely exponential. 
The important question however is whether the mode can be amplified significantly before its 
wavelength becomes small enough to be affected by ohmic resistivity. This, of course, depends 
on the magnetic Reynolds number. 

We have seen above that the so-called magnetorotational or Balbus-Hawley instability can 
develop in any disk in which the angular velocity decreases outward. In principle there is no 
other condition. However, it may be that the scale of the modes which are unstable according 
to (^9]) do not fit into the disk, i.e. they have a wavelength larger than the disk semithickness 
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H. This is the case if va/^ > H. Since in a thin disk VtH ~ Cg (see § p.3.1| ), the disk will be 
stable if va> Cs. 

Another condition which is implicit in the above presentation is, of course, that the magnetic 
field be coupled to the fluid. This may not be the case everywhere in protostellar disks, which 
are rather cold and dense (Gammie 1996; Fromang, Terquem & Balbus 2001). 

Ohmic resistivity can prevent the development of the instability if the scale on which it acts is 
comparable to the wavelength of the unstable modes. The time it takes for the magnetic field 
to diffuse over a scale ~ 1/k due to the effect of an Ohmic resistivity 77^ is ~ 1/(77^^^)- Since 
in the absence of resistivity it would grow on a timescale ~ l/f^, we see that Ohmic dissipation 
prevents the instability if tjb ~ Q/k"^, i.e. if the magnetic Reynolds number is of order unity. 
For more details, the reader is referred to Balbus & Blaes (1994), Jin (1996) and Papaloizou & 
Terquem (1997). Note that Fleming, Stone & Hawley (2000) have shown that even though the 
linear instability is not affected at higher values of the magnetic Reynolds number, the nonlinear 
instability cannot be sustained for values as high as lO'^ in some circumstances, which depend 
on the field topology for instance. However, the issue is not yet settled as other effects like 
Hall electromotive forces may resurrect the instability in regions where Ohmic resistivity acting 
alone would inhibit it (Wardle 1999, Balbus & Terquem 2001). 



4.3 Angular momentum transport in magnetized disks 

Numerical simulations have shown that this instability puts the energy it extracts from the 
disk differential rotation into fluctuations which transport angular momentum outward (see 
Balbus & Hawley 1998 and references therein). As noted above, transport of momentum is 
not a property of all types of turbulence. Furthermore, when it happens, it is not always 
against the velocity gradient, i.e. outward in an accretion disk. That has to be the case in an 
isolated dissipative system, but not if the energy driving the turbulence has an external source. 
We have already mentioned that the transport associated with thermal convection appears 
to be very weak compared to that associated with shear turbulence. In addition, there are 
strong indications that convection driven by external heating transports angular momentum 
inward in accretion disks (Ryu & Goodman 1992, Cabot 1996, Stone & Balbus 1996). In the 
case of the magnetorotational instability though, it can be shown that in the linear regime 
the transport is outward (Balbus & Hawley 1992a). This still holds in the nonlinear regime, 
as expected since the energy source for the turbulence is internal. Numerical simulations also 
show that most of the transport is done by the (magnetic) Maxwell stress, which dominates over 
the (hydrodynamic) Reynolds stress. Furthermore, magnetic fields are regenerated through a 
dynamo action so that the mechanism can sustain itself in an isolated system. 

It is important to realize that this instability and the resulting turbulence do not tend to make 
the disk angular velocity uniform. This is because the angular velocity is imposed by the 
gravitational potential of the object around which the disk rotates. The turbulence can smear 
the shear out only on scales smaller than the disk semithickness H. To smear the shear out 
on a scale /, the turbulent velocity Vt has indeed to be on the order of the shear over /, i.e. 
l\rdQ/dr\ ~ IQ, where Q is the angular velocity in the disk, li I ^ H, then Vt ~ c^, as in a thin 



disk Cs ~ HQ (see § |5.3.1| ). Therefore, only supersonic turbulence could modify the shear on 
scales larger than H. However, such a turbulence could not be maintained as the supersonic 
fluctuations would dissipate into shocks. The state toward which the disk evolves as a result 
of the turbulence is not a state of uniform rotation but one where all the mass is at the center 
and all the angular momentum at infinity (Lynden-Bell & Pringle 1974). 
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Long before a mechanism for producing turbulence in accretion disks had been identified, 
Shakura & Sunyaev (1973) proposed a prescription for modehng turbulent disks. We now 
describe this prescription, and discuss its validity in the context of magnetic turbulence. We 
then develop disk models based on this so-called a prescription. 



5 a disk models 

5.1 Evolution of turbulent disks 

Here we consider an axisymmetric disk rotating around a central object. We suppose the 
motion is in the plane of the disk, or, equivalently, we use the vertically averaged equations of 
mass conservation and motion. The velocity is v = {ur^r^ + «</>). The term rfi is the circular 
velocity around the central mass, and (ur, u^ are the components of the fluctuation velocity. 
Note that as the disk accretes, there is a net radial drift and the mean value of u.r is non zero. 
Here, the equation of mass conservation (j^) and the azimuthal component of the equation of 
motion (|) are: 

as 1 9 , , 

^ + -^(rE„.)=0. (54) 

where S is the surface mass density, i.e. the mass density integrated over the disk thickness. 
Here, according to the discussion in § f^.4|, we have neglected the viscous force arising from 



molecular viscosity. Multiplying equation (|55|) by r and using (|5^) , we obtain the angular 
momentum equation: 

^ (rE [r^ + u^)) + ^^ (r^S (rfi + u^) u,) = 0. (56) 

As pointed out by Balbus & Papaloizou (1999), to get a diffusion equation describing the disk 
evolution we need to smooth out the fluctuations over radius. To do so, we average the above 
equation over a scale large compared to that of the fluctuations, but small compared to the 
characteristic scale of the flow (i.e. r). Equation (^) then yields: 

— fSr^fi) + -Tr fSr^^ < Ur > +^r^ < UrUs >) = 0, (57) 

ot ^ ■' r or ^ ' 

where the brackets denote the radial average and we have neglected <u^> compared to rfi in 
the time derivative. This is justified because \ <u^> \<^ r^ and the systematic, evolutionary 
time derivative oi < u^> is limited. In the radial divergence term however, both <Ur> and 
< UrU^ > are second order, and therefore we retain all the terms. This equation is the same 
as that describing a viscous flow with v = {ur,rQ) and a stress tensor o"r</, = — S < WrW,/, > as 
given by equation (^). 

There are two contributions to the flux of angular momentum: the term Sr^f2 < -u^ > is the 
mean angular momentum advected through the disk by the velocity fluctuations (because of 



the accretion of mass), whereas the term Sr < UrU^ > represents the angular momentum 
fluctuations transported by the velocity fluctuations. 



Using equation (^) averaged over radius, we can rewrite equation (l57| ) under the form: 



rS < Ur >= — 



d_ 
dr 



d 



ir'^) k (^^' < -^-^ >) • 



Using equation (Q) again to eliminate < -u^ >, this leads to the diffusion equation: 




r^fi 



-' d_ 
dr 



Sr < UrU^ > 



(5^ 



In a steady state, equation ( p4|) gives rS < Ur >= constant. Then the accretion rate 



M = -27rrS < Ur > (59) 

is constant through the disk. Integration of the angular momentum equation { ^7\) then yields: 



^ 27r 



r 



(60) 



where Ri is the disk inner boundary. We have assumed here that the turbulent stress < UrUfp > 
vanishes at the disk inner edge (i.e. there is no torque at the boundary) and that the disk 
is Keplerian, so that f2 oc r~^'^. The above equation shows that for the mass to be accreted 
inward (i.e. < Ur > negative), the flux of angular momentum due to the fluctuations has to be 
positive, i.e. the fluctuations have to transport angular momentum outward. 

5.2 The a prescription 

We pointed out that the angular momentum equation ( ^Tl) is analogous to that describing a 
viscous flow with v = {ur,rQ) and a stress tensor ar^) = — S < UrU^ >. Therefore, it is 
tempting to push the analogy further and express the turbulent stress — S < UrU(f, > as if it 
derived from an enhanced 'turbulent viscosity' u, defined such that (equation [T^): 



In a Keplerian disk, this gives: 



L < UrU^ >= Lz/r— — . 
dr 



< UrU(f, >~ uQ. 



(61) 



(62) 



The equations presented in the previous section have been derived by Lynden-Bell & Pringle 
(1974; see also Pringle 1981) assuming a viscous flow with this expression of the stress tensor. 

The prescription proposed by Shakura & Sunyaev (1973) consists in writing z/ = vtH, where 
H is the disk thickness, assumed to be the maximum scale of the turbulent cells, and Vt is the 
turbulent velocity. They further define 

a = — , 

Co 
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where Cg is the sound velocity. Note that a < 1, otherwise the fluctuations would dissipate into 
shocks in such a way as to restore Vt < Cg. Equation ( |62|) can then be rewritten under the form: 



< UrUcj) >~ aCg. (63) 

Here we have used the fact that in a thin disk H ~ Cg/fi (see § [5.3.1| ). 

So far we have focussed on non magnetized disks. In these, there are strong indications that 
< UrUcf) > is either zero or negative. Magnetism is needed to correlate the velocities. The above 
discussion does apply to magnetized disks provided we replace < «,.«(/> > by < UrU^ — UArUA<j> >, 
where {uav, ua^) are the components of the fluctuations of the Alfven velocity (Shakura & 
Sunyaev 1973, Balbus & Hawley 1998). The extra term represents the Maxwell stress. 

The validity of the a prescription in the context of magnetic turbulence was discussed by Balbus 
& Papaloizou (1999). They first pointed out that, as long as < UrU^ > (or < UrU^ — UArUAip >) 
is positive, the disk dynamics is the same as if it were evolving under the action of a viscosity. 
In that case indeed, the diffusion coefficient in equation (|58|) is positive. We can then always 



define an a parameter according to (p5|), although a may not be constant through the disk. 

Note however that this implicitly assumes it is possible to average the equations over radius 
in the way described in § |5.1| . If the scale of the fluctuations and the characteristic disk scale 
are not well separated, such an average cannot be performed. Since the maximum scale of the 
fluctuations is of order the disk thickness if, this procedure requires that there is a scale large 
compared to H and small compared to r. This condition may be only marginally fulflUed in 
protostellar disks, in which H may be up to 0.1-0.2r. 

Balbus & Papaloizou (1999) further noted that for the a prescription to apply, the disk had to 
behave viscously not only in its dynamics but also in its energetics. The key point here is that a 
viscous disk dissipates locally the energy it extracts from the shear, whether in a steady state or 
not. This may not be the case in a turbulent disk where, if the turbulent cascade is not efficient, 
part of the energy may be advected with the flow. As we have not addressed the energetics of 
viscous disks above, we will not go into the details of the discussion here. These can be found in 
Balbus & Papaloizou (1999), who showed that in disks subject to MHD turbulence the energy 
extracted from the shear is indeed dissipated locally (through the turbulent cascade) whether 
the disk is evolving or not. Note that this is in general not the case when the turbulence is due 
to self-gravitating instabilities. In that case indeed, part of the energy is transported by waves 
through the disk. 

Most theoretical protostellar disk models have relied on the a parametrization. We have com- 
mented that MHD turbulence does lend itself to this prescription in thin disks. However, 
because MHD instabilities develop only in an adequately ionized fluid, they may not operate 
everywhere in protostellar disks (Gammie 1996). Therefore it is likely that the parameter a is 
not constant through these disks. It may even be that only parts of these disks can be described 
using this a prescription. However, we may still learn about disks from these models in the 
same way as we learned about stars from simple polytropic models. Therefore, we now describe 
these models. 



5.3 Vertical structure 

In a thin disk, the radial gradients of temperature and pressure are much smaller than the 
vertical gradients. The vertical and radial structures therefore decouple, and can be calculated 
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separately. Here we present the calculation of the vertical structure, i.e. we determine pressure, 
temperature and density as a function of height z. Once this is done, given an initial surface 
density distribution, the radial structure can be computed using the diffusion equation (0). 
The section below is based on Papaloizou & Terquem (1999), which the reader is referred to 
for more details. 



5.3.1 Basic equations 

The equations describing the disk vertical structure are the equation of vertical hydrostatic 
equilibrium: 

i£^._^i^, (64) 

pdz (^2 + ^2)3/2' 

and the energy equation, which states that the rate of energy removal by radiation is locally 
balanced by the rate of energy production by viscous dissipation: 

— = -pun , (65) 

where JF is the radiative flux of energy through a surface of constant z which is given by: 

-IGaT^dT 
Sup dz ' 

Here G is the gravitational constant, M^, is the central mass, P is the pressure, T is the 
temperature, k is the opacity, which in general depends on both p and T, and a is the Stefan- 
Boltzmann constant. 

To close the system of equations, we relate P, p and T through the equation of state of an ideal 
gas: 

P=^, (67) 

prriH 

where k is the Boltzmann constant, p is the mean molecular weight and rriH is the mass of 
the hydrogen atom. If we limit our calculations to temperatures lower than 4,000 K (typical of 
protostellar disks), then, at the densities of interest, hydrogen is in molecular form. Since the 
main component of protostellar disks is hydrogen, we then have p = 2. We note that for the 
temperatures we consider transport of energy by convection can be neglected. 

We denote the isothermal sound speed by Cg {cj. = P/ p). We adopt the a-parametrization of 
Shakura & Sunyaev (1973), so that the kinematic viscosity is written v = ac^/Q (equations ^^ 
and ||6^). In general, a is a function of both r and z. With this formalism, equation ( |65D 
becomes: 

dT 9 

- ^ -„0P. (68) 

Note that, in a thin disk, z <^ r. Then, if we set \dP/dz\ ~ P/H and z ~ if in equation (|6^), 
we get Cg ~ QH, where Q is the Keplerian velocity [Q"^ = GM^/r^) and we have used c^ = P/p. 
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The thin disk approximation then requires c^ <^ rQ, i.e. the angular velocity to be highly 
supersonic. 

5.3.2 Boundary conditions 

We have to solve three first order ordinary differential equations for the three variables JF, P 
(or equivalently p), and T as a function of z at a given radius r. Accordingly, we need three 
boundary conditions at each r. We shall denote with a subscript s values at the disk surface. 

The first boundary condition is obtained by integrating equation (^) over z between —H and 
H. Since by symmetry JF(2; = 0) = 0, this gives: 

J^s = ^Mstn\ (69) 

OTT 

where we have defined Mgt = 37r(z/)S, with < u > being the vertically averaged viscosity. If 
the disk were in a steady state, Mgt would not vary with r and would be the constant accretion 
rate through the disk (eq. pD| and p2[ with Ri <^ r). In general however, this quantity does 



depend on r. 

Another boundary condition is obtained by integrating equation (|^) over z between H and 
infinity. A detailed derivation of this condition is presented in the Appendix A of Papaloizou 
& Terquem (1999). Here we just give the result: 

Ps = -, (70) 

Kg 

where Tab is the optical depth above the disk. This condition is familiar in stellar structure, 
where Q^H would be replaced by the acceleration of gravity at the stellar surface. Since we 
have defined the disk surface such that the atmosphere above the disk is isothermal, we have 
to take Tab ^ 1- Providing this is satisfied, the results do not depend on the value of Tab we 
choose. 

A third and final boundary condition is given by the expression of the surface temperature (see 
Appendix A of Papaloizou & Terquem 1999 for a detailed derivation of this expression): 

2a {t^ - Tf) - '— - -Mg,Q' = 0. 71 

Here the disk is assumed immersed in a medium with background temperature T^. The surface 
opacity k^ in general depends on both Tg and ps and we have used c^ = kT /{pniH)- The 
boundary condition (^) is the same as that used by Levermore & Pomraning (1981) in the 
Eddington approximation (their eq. [56] with 7 = 1/2). In the simple case when T^ = and the 
surface dissipation term involving a is set to zero, with Mgt being retained, it simply relates 
the disk surface temperature to the emergent radiation flux. 

5.3.3 Disk models 

At a given radius r and for a given value of the parameters M^t and a (which uniquely determine 
the disk model), the disk vertical structure is obtained by solving equations (0), (|66|) and ( |68|) 
with the boundary conditions (p9|), ( |70| ) and (^TJ). 
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Numerical results are published, e.g., by Papaloizou & Terquem (1999) (see also Lin & Pa- 
paloizou 1980, 1985; Bell & Lin 1994; Bell et al. 1997), in which the opacity is taken from 
Bell & Lin (1994). This has contributions from dust grains, molecules, atoms and ions. It is 
written in the form k = K,ip°'T'' where /tj, a and b vary with temperature. For a specified M, 



stt 



they calculate the value of H, the vertical height of the disk surface, iteratively. Starting from 
an estimated value of H, after satisfying the surface boundary conditions, the equations are 
integrated down to the mid-plane z = Q. The condition that JF = at 2; = is not in general 
satisfied. An iteration procedure is then used to adjust value of H until jF = 0at2; = 0toa 
specified accuracy. 
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